work.dir = "/homes/26/jc184963/climate_change2/spp_richness/all.spp.full/"; setwd(work.dir)
km.dir = "/homes/26/jc184963/climate_change2/spp_richness/all.spp.null/"
library(SDMTools)
out.dir = "/homes/26/jc184963/climate_change2/spp_richness/change.figs/full.dist.w.clipped.1990/"

asc.A1B.1990 = read.asc.gz(paste(km.dir,"A1B.ALL.SPP.no.disp.1990.asc.gz",sep=''))
asc.A1B.2080 = read.asc.gz("A1B.ALL.SPP.full.disp.2080.asc.gz")

asc.A2.1990 = read.asc.gz(paste(km.dir,"A2.ALL.SPP.3km.disp.1990.asc.gz",sep=''))
asc.A2.2080 = read.asc.gz("A2.ALL.SPP.full.disp.2080.asc.gz")

asc.B1.1990 = read.asc.gz(paste(km.dir,"A1B.ALL.SPP.3km.disp.1990.asc.gz",sep=''))
asc.B1.2080 = read.asc.gz("B1.ALL.SPP.full.disp.2080.asc.gz")

change.A1B = asc.A1B.2080 - asc.A1B.1990
change.A2 = asc.A2.2080 - asc.A2.1990
change.B1 = asc.B1.2080 - asc.B1.1990

write.asc.gz(change.A1B,paste(out.dir,"change.A1B.null.original.dist.asc",sep=""))
write.asc.gz(change.A2,paste(out.dir,"change.A2.asc",sep=""))
write.asc.gz(change.B1,paste(out.dir,"change.B1.asc",sep=""))

output.dir ="/homes/26/jc184963/climate_change2/spp_richness/change.spp.rich.image/full.dist.clipped.1990/"
change.dir = out.dir; setwd(change.dir)
 tasci=list.files()
 
pnts=cbind(x=c(148,149.3,149.3,148), y=c(-13.2,-13.2,-17.2,-17.2))
bins = seq(0,1,length=112)
#cols = c(rep('darkolivegreen3',bins),colorRampPalette(c("yellow","orange","red"))(100)[bins:100])
cols = c(colorRampPalette(c("red","orange","yellow","green","turquoise","blue","purple"))(100))
cols = c(colorRampPalette(c("red","orange","yellow","green","blue","purple","black"))(100))
tasc=read.asc.gz("change.A1B.asc.gz")
        png(paste(output.dir,"A1B.spp rich change.png", sep=""),bg= "white",height=800, width=500)
        image(tasc, col=cols,xlab="", ylab="",font.main=3)            #,zlim=c(0,1) , main= "spp rich change A1B"
        legend.gradient(pnts,cols=c(cols),limits=c(-62,50), title='Spp richness change')
        dev.off()

tasc2=read.asc.gz("change.A2.asc.gz")
        png(paste(output.dir,"A2.spp rich change.png", sep=""),bg= "white",height=800, width=500)
        image(tasc2, col=cols,xlab="", ylab="",font.main=3)            #,zlim=c(0,1) , main= "spp rich change A1B"
        legend.gradient(pnts,cols=c(cols),limits=c(-62,50), title='Spp richness change')
        dev.off()

tasc3=read.asc.gz("change.B1.asc.gz")
        png(paste(output.dir,"B1.spp rich change.png", sep=""),bg= "white",height=800, width=500)
        image(tasc3, col=cols,xlab="", ylab="",font.main=3)            #,zlim=c(0,1) , main= "spp rich change A1B"
        legend.gradient(pnts,cols=c(cols),limits=c(-62,50), title='Spp richness change')
        dev.off()

  zlimits = range(c(tasc[,],tasc2[,],tasc3),na.rm=TRUE)

       pnts=cbind(x=c(147.8,149.3,149.3,147.8), y=c(-13.2,-13.2,-17.6,-17.6))
png(paste(output.dir,"spp rich change.full.png", sep=""),width=17,height=5,units='cm',res=600,pointsize=6,bg= "white")
 	 par(mfrow = c(1,3),oma=c(3,3,0.5,0.5),mar=c(0,0,0,0),cex=.8,cex.axis=1.5)
     image(tasc3, axes=FALSE,ann=FALSE,col=cols,zlim=c(-44,128))
     image(tasc, axes=FALSE,ann=FALSE,col=cols,zlim=c(-44,128))
     image(tasc2, axes=FALSE,ann=FALSE,col=cols,zlim=c(-44,128))
     legend.gradient(pnts,cols=c(cols),limits=c(-44,128), title='Rich change')
     mtext("B1",side=1,at=.14,cex=1.5,line=1.5,outer=TRUE)
          mtext("A1B",side=1,at=.49,cex=1.5,line=1.5,outer=TRUE)
               mtext("A2",side=1,at=.81,cex=1.5,line=1.5,outer=TRUE)
    dev.off()


###############species richness fancy-like #####################################
work.dir = "/homes/26/jc184963/climate_change2/spp_richness/change.figs/full.dist.w.clipped.1990/"; setwd(work.dir)
output.dir ="/homes/26/jc184963/climate_change2/spp_richness/change.spp.rich.image/full.dist.clipped.1990/"
library(SDMTools)
 tasci=list.files()

tasc=read.asc.gz("change2.A1B.asc.gz")
tasc2=read.asc.gz("change2.A2.asc.gz")
tasc3=read.asc.gz("change2.B1.asc.gz")

        pnts=cbind(x=c(147.8,149.3,149.3,147.8), y=c(-13.2,-13.2,-17.6,-17.6))
 tt = tasc[which(is.finite(tasc))];tt=as.numeric(as.character(cut(tt,breaks=seq(-44,128,4),labels = seq(-43.5,127.5,4)))) ;tasc[which(is.finite(tasc))] = tt
 tt = tasc2[which(is.finite(tasc2))];tt=as.numeric(as.character(cut(tt,breaks=seq(-44,128,4),labels = seq(-43.5,127.5,4)))) ;tasc2[which(is.finite(tasc2))] = tt
 tt = tasc3[which(is.finite(tasc3))];tt=as.numeric(as.character(cut(tt,breaks=seq(-44,128,4),labels = seq(-43.5,127.5,4)))) ;tasc3[which(is.finite(tasc3))] = tt
 cols = c(colorRampPalette(c('red','yellow'))(5),'grey',colorRampPalette(c('lightblue','darkblue'))(15))
png(paste(output.dir,"spp rich change.full.stdzd.png", sep=""),width=17,height=5,units='cm',res=600,pointsize=6,bg= "white")
 	 par(mfrow = c(1,3),oma=c(3,3,0.5,0.5),mar=c(0,0,0,0),cex=.8,cex.axis=1.5)
     image(tasc3, axes=FALSE,ann=FALSE,col=cols,zlim=c(-44,128))
     image(tasc, axes=FALSE,ann=FALSE,col=cols,zlim=c(-44,128))
     image(tasc2, axes=FALSE,ann=FALSE,col=cols,zlim=c(-44,128))
     legend.gradient(pnts,cols=c(cols),limits=c(-44,128), title='Rich change')
     mtext("B1",side=1,at=.14,cex=1.5,line=1.5,outer=TRUE)
          mtext("A1B",side=1,at=.49,cex=1.5,line=1.5,outer=TRUE)
               mtext("A2",side=1,at=.81,cex=1.5,line=1.5,outer=TRUE)
    dev.off()

